ST310 Course Project
Chris Chia, Mun Fai Chan, Zhen Yen Chan
29/04/21
ST310 Course Project
The aim of this project is to demonstrate the use of machine learning techniques learnt in the ST310 Machine Learning module. We rely on both linear (Logistic Regression) and non-linear methods (Random Forests and Gradient Boosted Decision Trees) for a classification problem.
Remark: If executing this notebooks as a .Rmd file, ensure that all libraries and dependencies are installed and that the chunks are executed in sequential order.
Dataset
We have obtained data from the Kaggle March Tabular Playground competition. The data consists of anonymised features, which correspond to a binary outcome variable; in other words the task is a classification problem. Although the data is anonymised, the challenge description states:
The dataset is used for this competition is synthetic but based on a real dataset and generated using a CTGAN. The original dataset deals with predicting the amount of an insurance claim. Although the features are anonymized, they have properties relating to real-world features.
The dataset consists of 30 features, 19 categorical variables and 11 numerical variables., and a binary outcome variable target. 73.5% of the observations have ‘target==1,’ suggesting an imbalanced dataset. We will need to process these categorical variables in some manner, which we shall consider in subsequent sections.
library(tidyverse)
library(ggplot2)
library(patchwork)
library(broom)
library(tidymodels)
library(arm)
library(car) #outlier
library(tidyr)
library(glmnet)
library(caret)
library(randomForest)
library(xgboost)
library(catboost)
library(e1071)
library(data.table)
library(doParallel)
library(foreach)
library(pROC)
library(pdp)
df <- read.csv(file = "../data/train.csv")
cat(c("The dimensions of the array are :", dim(df)[1], ", ", dim(df)[2]))The dataset consists of 30 features, 19 categorical variables and 11 numerical variables., and a binary outcome variable target. We will need to process these categorical variables in some manner, which we shall consider in the subsequent section.
A challenge is that we have no specific domain knowledge about the meaning of particular features, and whether they may be useful. Nonetheless, we can resort to model interpretability methods (Molnar, Christoph 2019) to obtain an ex-post explanation of our models. In the linear case, we can examine coefficients and their statistical significance, and in the case of tree-based models, we can utilise feature importances, partial dependence plots, and Shapley Values.
Preprocessing
We first remove the first column id, which represents an unique identifier for each observation. We then convert the first 19 columns, which are categorical variables into the factor data type in R, To reduce computational time, we subsample 4000 observations from the complete dataset of 300,000 observations. For our subsequent analysis, we will use only this subsample of 4000 observations, although the approaches can be extended to a larger dataset given greater computational resources. We further partition the 4000 observations, into a train set of 3000 observations and a test set of 1000 observations. All training and tuning are performed on the training data, and we will consider model performance on the test set performance scores.
Remark: Different models would require different preprocessing for the categorical variables, and thus we will further address in subsequent sections; for example the glm package can accept factors as a data type, whereas for glmnet they must be converted to the model.matrix format.
Remark: When encoding the categorical variables as dummy variables in our subsequent analysis, we find that there are over 600 unique categories across all 19 of the categorical variables. We could also consider different methods to process categorical variables, such as merging less frequently occuring categories, or target encoding (Parr, Terence and Howard, Jeremy 2019, ch. 6). Nonetheless, we stick with the factors and dummies approach given our lack of familiarity with these other methods.
cat_feats = 1:19 # the first 19 columns are categorical
cont_feats <- 20:30 # and the next 11 are continuous
target_col <- 31 # target
# convert cats to factors
df <- column_to_rownames(df, var = "id") %>%
mutate_if(is.character,as.factor) %>%
mutate_at(vars(target), factor)
# subsample the data for faster model imputation
set.seed(1)
sam = sample(1:nrow(df), 4000)
df_sample = df[sam,]
# Partition data into train and test; test will be our oos data
set.seed(1)
df_split <- initial_split(df_sample, prop = 3/4)
df_train <- training(df_split)
df_test <- testing(df_split)Exploratory Data Analysis
Univariate EDA
We observe that the univariate distributions of the continuous variables are all multi-modal and non-normal, but they are all normalised to the range of \([0, 1]\).
# flatten df into using pivot_longer and plot distribution
df %>% pivot_longer(cols = starts_with("cont"), names_to = "cont") %>%
ggplot(aes(x = value))+
geom_histogram(bins = 100, alpha = 0.85)+
ggtitle("Continuous features distribution")+
facet_wrap(cont~.,scales = "free") +
theme_minimal()
From the distributions of categorical variables, we see that there are variables with substantially more observations in one category, and also variables with a high number of (\(>50\)) categories, which will be an issue to address in our preprocessing step.
# flatten df into using pivot_longer and plot distribution
df %>% pivot_longer(cols = contains(c("cat", "target")), names_to = "cat") %>%
ggplot(aes(x = value))+
geom_bar(alpha = 0.85)+
ggtitle("Categorical features distribution")+
facet_wrap(cat~.,scales = "free", ncol = 4) +
theme_minimal(base_size = 30)
Bivariate EDA
We group the continuous variables by the target and plot them as boxplots to check for any obvious differences discernible by eye. From the plots, claims with target == 1 have lower values of cont3 on average (median) than claims with target == 0, hence we would expect a negative relationship between cont3 and target. Claims with target=1 also have a higher median value of cont4 than claims with target == 0.
# flatten df into using pivot_longer
# group by target and plot distribution
df[, c(cont_feats, target_col)] %>%
pivot_longer(cols = starts_with("cont"), names_to = "var", values_to="value") %>%
ggplot(aes(x=target,y=value), fill=factor(value)) +
geom_boxplot() + coord_flip() + facet_wrap(~var, scales="free_x")
From our conditional boxplots, we can identify potential outliers. In particular we can spot many potential outliers for cont0, cont5, cont7, cont8, cont9 and cont10. To investigate this, we identified observations that lie at the extreme percentiles of these variables. Using the Hampel filter, which considers points lying outside the median plus or minus 3 mean absolute deviations as outliers, more than 200 observations per variable were classified as such. This suggests that these may not be outliers but that the distribution is just heavy-tailed. Without further information on the reasonable scale of values that individual variables can take (along with the fact that they are all normalised), we find it challenging to clearly identify outliers through descriptive statistics and decide not to exclude any such points using this approach. We will proceed to detect for outliers using a model approach in a subsequent section.
hampel_filter <- function(df){
lower_bound <- median(df) - 3 * mad(df, constant = 1)
upper_bound <- median(df) + 3 * mad(df, constant = 1)
outlier_ind <- which(df < lower_bound | df > upper_bound)
return(outlier_ind)
}
percentile_filter <- function(df, lq = 0.001, uq = 0.999){
lower_bound <- quantile(df, lq)
upper_bound <- quantile(df, uq)
outlier_ind <- which(df < lower_bound | df > upper_bound)
return(outlier_ind)
}
hampel_count <- function(x){length(hampel_filter(x))}
pct_count <- function(x){length(percentile_filter(x))}
outlier_counts <- df_train[, cont_feats] %>% map_dfr(hampel_count)
outlier_counts[2, ] <- df_train[, cont_feats] %>% map_dfr(pct_count)
outlier_countsFor categorical variables, we use stacked bar plots to show the percentages of observations in each category that correspond to target == 0 and target == 1 respectively. A much larger proportion of claims with cat13 == B appear to be associated with target == 1 compared to cat13 == A. On the other hand, a much larger percentage of claims correspond to target == 0 if cat18 is A or B, than if cat18 is C or D.
# flatten df into using pivot_longer
# group by target and plot distribution
df[, c(cat_feats, target_col)] %>%
pivot_longer(cols = starts_with("cat"), names_to = "cat", values_to="value") %>%
ggplot(aes(x = value, fill=target)) +
geom_bar(position="fill") +
scale_y_continuous(name = "Within group Percentage", labels = scales::percent) +
facet_wrap(~cat, scales="free_x", ncol = 4) +
theme_minimal(base_size = 40)
We also inspect the correlation matrix for our continuous variables. There seems to be a cluster of variables - cont1, cont2, cont8 - that are highly correlated with each other. This could potentially lead to problems with multicollinearity when using linear models. In addition, we also consider the pearson correlation and our continuous variables, although not necessarily meaningful in this case given our target is binary, and in addition given the multi-modal distributions of the continuous variables we identified in our univariate exploratory data analysis. However, we note that the continuous variables have pearson correlations of roughly \([-0.2, 0.2]\), suggesting that there is some signal between the features and the target. Finally, we note that the pearson correlation only describes a linear and pairwise association between the variables; as such there could be complex non-linear associations and interactions between the variables as well.
cor_matrix <- cor(df[, cont_feats])
heatmap(cor_matrix, main="Correlation Matrix (Clustered)")
cor_with_target <- rownames_to_column(data.frame(cor(df[, c(cont_feats)], as.numeric(df$target))))
names(cor_with_target) <- c("feature", "pearson_corr")
cor_with_target %>%
ggplot(aes(x=reorder(feature, -pearson_corr), y = pearson_corr)) +
geom_bar(stat='identity') + coord_flip() + ggtitle("Pearson Correlation of Continuous with Target")
We can use Principal Components Analysis (PCA) as a means to visualise the data in low-dimension, to determine if there are any explicitly discernible trends. By eye, the classes do not appear to be linearly separable - which suggest a non-linear method may be more effective. There do not appear to be any significant difference in the PCA representations for each class, although there are many observations with target == 1 closer to the bottom of the ellipsoid formed by the first two PCA loadings. We note that there is a connection between PCA and Ridge (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome 2009), and also that the first two principal components only capture about \(\approx 60\%\) variation in the continuous variables.
pcs <- prcomp(df[,cont_feats])
set.seed(2021)
data.frame(pc1=pcs$x[,1], pc2=pcs$x[,2], target=df[, "target"]) %>%
ggplot(aes(x = pc1, y = pc2, colour = target)) +
geom_jitter(alpha=0.7) + ggtitle('Principal Components')
# cumulative variance
cumul_var <- cumsum(pcs$sdev^2 / sum(pcs$sdev^2))
ggplot(data.frame(feature = 1:11, cumul_var = cumul_var)) +
geom_line(aes(x = feature,y = cumul_var)) + ggtitle("Cumulative Explained Variance Ratio") +
scale_y_continuous(labels=percent)
Modelling
We first consider a naive model that always predicts a single label (in this case \(0\)) given it is the most frequently occurring class. In this case, the accuracy would \(1 - \hat{y} = 0.745\) This illustrates the issue with the accuracy metric - the naive accuracy is high due to the nature of the data, and hence subsequent model performances need to be compared with this benchmark.
We also consider the ROC-AUC metric (Receiving Operator Characteristic Area Under the Curve), and accuracy metric. Given the problem is one related to insurance, the modelling outcome of interest is not only to obtain the correct predictions (accuracy), but also accurate probabilities, which is what the AUC metric measures. AUC measures the performance of the classifier across all possible decision thresholds, which makes it more helpful to compare performance than the accuracy metric with default decision threshold set to 0.5. As we do not know the specific threshold relevant to the insurance context of the problem, we will report both accuracy and ROC for all models, and compare it to that of the baseline model in [section 1]. On the other hand, predicting all 0s would result in the lowest ROC-AUC score of 0.5 as expected.
# train-test
X_train = as.matrix(df_train[, cont_feats])
y_train = as.numeric(as.matrix(df_train$target))
X_test = as.matrix(df_test[, cont_feats])
y_test = as.numeric(as.matrix(df_test$target))
diagnosis <- function(train_pred, test_pred, train_true, test_true){
train_classes <- ifelse(train_pred > 0.5, 1,0)
test_classes <- ifelse(test_pred > 0.5, 1,0)
acc1 <- mean(train_classes == train_true)
auc1 <- auc(roc(train_true, train_pred, quiet=TRUE))
acc2 <- mean(test_classes == test_true)
auc2 <- auc(roc(test_true, test_pred, quiet=TRUE))
data.frame(train_acc=acc1, train_auc=auc1, test_acc = acc2, test_auc=auc2)
}
results = data.frame(diagnosis(rep(0, dim(df_train)[1]),
rep(0, dim(df_test)[1]),
df_train$target, df_test$target),
row.names=c("naive"))
resultsLogistic Regression (SGD)
To fulfill the project requirements, we demonstrate a `from scratch’ Stochastic Gradient Descent routine for Logistic Regression. The deriviation follows (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome 2009, 120–26)
The logistic loss is given by: \[l(\boldsymbol{\beta}) = -\sum_{i = 1}^{N} y_{i} log(p(x_{i} ; \boldsymbol{\beta})) + (1 - y_{i}) log(1 - p(x_{i} ; \boldsymbol{\beta})) = \sum_{i = 1}^{N} \left [ y_{i}log \left (\frac{p(x_{i} ; \boldsymbol{\beta})}{1 - p(x_{i} ; \boldsymbol{\beta})} \right) + log(1 - p(x_{i} ; \boldsymbol{\beta})) \right ]\]
\[ l(\boldsymbol{\beta})= -\sum_{i = 1}^{N}\left [y_{i} \boldsymbol{\beta}^{T} x_{i} - log(1 + exp(\boldsymbol{\beta}^{T}x_{i})) \right ]\]
With the inclusion of a reularisation term, in this case a \(L^{2}\) penalty:
\[ l(\boldsymbol{\beta}) = -\sum_{i = 1}^{N} \left [y_{i} \boldsymbol{\beta}^{T} x_{i} - log(1 + exp(\boldsymbol{\beta}^{T}x_{i})) \right ] - \lambda \boldsymbol{\beta}^{T} \boldsymbol{\beta}\]
The gradient of the loss function is given by:
\[\nabla(\boldsymbol{\beta}) = -\sum_{i = 1}^{N} \left [ y_{i} x_{i} - \frac{x_{i}exp(\boldsymbol{\beta}^{T} x_{i})}{1 + exp(\boldsymbol{\beta}^{T} x_{i})} \right] - \lambda 2\boldsymbol{\beta}\]
We use the Barzilai-Borwein method (Murphy, Kevin P. 2012, 444–45) to determine the step size.
# binary crossentropy / log-loss
log_loss <- function(x, y, betas, lambda){
logits <- x %*% betas
- (t(y) %*% logits - sum(log(1 + exp(logits))) + lambda * t(betas) %*% betas) / dim(x)[1]
}
# logistic regression gradients
gradients <- function(x, y, betas, lambda){
logits <- x %*% betas
- (t(x) %*% (y - exp(logits)/(1 + exp(logits)))) - lambda *2 * betas / dim(x)[1]
}
p = dim(X_train)[2]
lambda = 0
n_iters <- 100
init_step_size <- 1e-6
set.seed(2021)
beta_init <- matrix(rnorm(p),nrow=p)
beta_path <- matrix(rep(0, n_iters * p), nrow = n_iters, ncol=p)
beta_path[1,] = beta_init
last_grad <- grad <- gradients(X_train, y_train, beta_path[1,], lambda)
beta_path[2,] = beta_init - init_step_size * grad
grad <- gradients(X_train, y_train, beta_path[2,], lambda)
losses <- rep(0, n_iters)
for (i in 3:n_iters){
step_size <- as.numeric(t(beta_path[i - 1,] - beta_path[i - 2,]) %*% (grad - last_grad) /
(t(grad - last_grad) %*% (grad - last_grad)))
beta_path[i,] <- beta_path[i - 1,] - step_size * grad
last_grad <- grad
grad <- gradients(X_train, y_train, beta_path[i, ], lambda)
losses[i] <- log_loss(X_train, y_train, beta_path[i,], lambda)
}
ggplot(data.frame(step = 3:n_iters, loss=losses[3:n_iters])) +
geom_line(aes(x = step, y = loss)) +
ggtitle("Binary Crossentropy vs. Iterations")
pred_train <- as.numeric(1 / (1 + exp(-X_train %*% beta_path[100,])))
pred_test <- as.numeric(1 / (1 + exp(-X_test %*% beta_path[100,])))
results["sgd",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["sgd",]Logistic Regression (Few Predictors)
As a baseline model, we build a simple logistic regression with a few predictors, which are selected from our exploratory data analysis to have discernible differences in target. These are cont3, cont4, cat13, cat18. We use the glm package to fit a logistic regression.
glm1 <- glm(target~cont3+cont4+cat13+cat18,data=df_train, family=binomial(link="logit"))
summary(glm1)
Call:
glm(formula = target ~ cont3 + cont4 + cat13 + cat18, family = binomial(link = "logit"),
data = df_train)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5707 -0.6702 -0.6021 0.3931 2.0095
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.4821 1.0756 -1.378 0.16824
cont3 -0.6554 0.2276 -2.880 0.00398 **
cont4 -0.4772 0.2063 -2.313 0.02070 *
cat13B 1.8681 0.3346 5.584 2.36e-08 ***
cat18B 0.5581 1.0714 0.521 0.60243
cat18C 2.6081 1.0816 2.411 0.01590 *
cat18D 2.4342 1.0805 2.253 0.02427 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3406.6 on 2999 degrees of freedom
Residual deviance: 2962.8 on 2993 degrees of freedom
AIC: 2976.8
Number of Fisher Scoring iterations: 4All four predictors are significant at 5% level, with the most significant predictors being… The coefficients of the dummy variables indicate the average difference between the log odds of that factor level group compared to the baseline level group. In our regression, the first category (A) of each categorical variable is taken as the baseline group. For example, the coefficient of cat13b implies the following equation: \[\text{logit} P(\text{target}=1|\text{cat13b} = 1) – \text{logit} P(\text{target}=0|\text{cat13b} = 1) = 1.8681. \] In words, this means that claims with cat13b = 1 have 33% higher log odds of target=1 than claims with \text{cat13b} == 0. However, we can only assert that there is some association - in this case it is statistically significant relationship with the p-value being 2.36e-08 - between the presence of the variable cat13b = 1 and the target being 1, but not whether there is a causation; in particular, given the lack of knowledge of what the variable represents, we cannot form any meaningful hypothesis.
To interpret the coefficients of the continuous variables, risk ratios need to be calculated using specific pairs of values of the predictors. The non-linearity of the logistic function means that a change of one unit in the value of a predictor is not the same across the range of the predictor.
pred_train <- predict(glm1, df_train, type="response")
pred_test <- predict(glm1, df_test, type="response")
results["glm-small",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["glm-small",]Logistic Regression (All Predictors)
We now run a logistic regression using all variables, using the glm package.
glm2 <- glm(target~., data=df_train, family=binomial(link="logit"),
control = list(maxit = 100))
# display(glm2)
pred_train <- predict(glm2, df_train, type="response")
glm2$xlevels = lapply(df[,cat_feats], levels)
pred_test <- predict(glm2, df_test, type="response")
results["glm-full",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["glm-full",]# anova(glm1, glm2, test="Chisq")Outlier Detection
We attempt to detect outliers under a logistic regression. Formally, outliers are defined as observations with a response vector that is unusual conditional on covariates. They are formally identified through studentized residuals. Intuitively, outliers have large residuals and we can formally test (by looking at the Bonferroni–adjusted p-values) if these residuals are significantly larger than the other observations. We use the outlierTest function from the car to determine the outliers from the glm model.
Observations that are far from the average covariate pattern are considered to have high leverage and can be measured using the hat value. Here, there are many points with high leverage.
outlierTest(glm2)
outliers <- as.numeric(names(outlierTest(glm2)$p))Finally, we measure for influence, which is an observation that is an outlier and have high leverage. These are likely to influence the regression coefficients and influence can be thought of as the product of leverage and outlier. Here, we plot studentised residuals against hat-values with the size of a circle being proportional to the Cook’s distance of an observation- a measure of influence.
influenceIndexPlot(glm2, vars = "hat")influencePlot(glm2)Here, we observe that there are a number of observations with high influence - outliers with high leverage. Thus, we remove these observations and compare the performance of our updated model with the original model.
In the later models, we will also exclude the same outliers.
influencers <- as.numeric(rownames(influencePlot(glm2)))
glm2_influencers <- update(glm2, subset = c(-influencers))
glm2_outliers <- update(glm2, subset = c(-outliers))
removal_list <- union(outliers, influencers)
glm2_removed <- update(glm2, subset = c(-removal_list))
compareCoefs(glm2, glm2_influencers, glm2_outliers, glm2_removed)
# actually just use glm2 and glm2_removedRegularised Logistic Regression
Given the issue of high dimensionality, we consider a regularised form of logistic regression. We use the glmnet package; in doing so we need to convert the data type into matrices.The glmnet package requires the data to be in a matrix data type, and hence we make the corresponding adjustment.
X_train = df_train[, -length(df_train)]
y_train <- df_train$target
X_test = df_test[, -length(df_test)]
y_test <- df_test$target
X_train = model.matrix(~., X_train)
X_test = model.matrix(~., X_test)
glm3 <- cv.glmnet(X_train, y_train,
family="binomial"(link="logit"), alpha=0)
glm3
Call: cv.glmnet(x = X_train, y = y_train, family = binomial(link = "logit"), alpha = 0)
Measure: GLM Deviance
Lambda Index Measure SE Nonzero
min 0.1422 80 0.7846 0.01822 458
1se 0.3606 70 0.8010 0.01636 458pred_train <- as.numeric(predict(glm3, X_train, type="response"))
pred_test <- as.numeric(predict(glm3, X_test, type="response"))
results["glm-ridge",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["glm-ridge",]Random Forest
To improve performance, we draw on the usage of non-linear tree-based models, specifically Random Forests (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome 2009). Intuitively, a random forest averages different decision trees (known as bagging) so as to reduce the variance of individual trees. We use the
rf_recipe <- recipe(target~., data = df_train)
rf_model <-
rand_forest(trees=100) %>%
set_engine("ranger", importance="impurity", seed=2021) %>%
set_mode("classification")
rf_workflow <- workflow() %>%
add_recipe(rf_recipe) %>%
add_model(rf_model)
rf1 <- fit(rf_workflow, df_train)
ranger_obj <- pull_workflow_fit(rf1)$fit
# plot feature importances
rf_feat_imp <- rownames_to_column(data.frame(ranger_obj$variable.importance));
names(rf_feat_imp) <- c("feat", "importance")
ggplot(rf_feat_imp, aes(x = reorder(feat, importance), y = importance))+
geom_bar(stat="identity", position="dodge")+ coord_flip()+
ylab("Feature Importance (Gini Impurity)")+
xlab("Feature")+
ggtitle("Random Forest Feature Importances")
#evaluate metrics
pred_train <- unlist(predict(rf1,df_train,type='prob')[,2])
pred_test <- unlist(predict(rf1, df_test, type='prob')[,2])
results["rf",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["rf",]temp <- function(x){partial(ranger_obj, train=df_train, pred.var = x, plot = TRUE, plot.engine = "ggplot2", paropts = list(.packages = "ranger"))}
plots <- lapply(names(df_train)[cont_feats], temp)
wrap_plots(plots)
The random forest has a test accuracy of 0.836 and a test AUC of 0.872, which is higher than the previous linear models. However, the train accuracy and AUC are significantly higher than the test performance, which suggest that there may be some degree of overfitting to the train data.
Random forests can be used to rank the importance of different features. Specifically, the x-axis is the Mean Decrease Accuracy, which reports how much accuracy the model loses when we exclude this variable. The more the accuracy falls by, the more important the particular variable is. Note here that for categorical variables, each level of the category is classified as a single variable. In this plot, we recorded the 30 most important variables.
We then run another random forest with the most important features. In doing so, we hope to reduce the degree of overfitting by reducing the complexity of the model. However, it does not make sense to drop some levels of a categorical variable while including the other levels. Hence, as long as a level is present in the top 30 features, we will include the entire category in our updated random forest model. This results in us keeping only 22 variables, from an initial 30.
Our reduced random forest has a slightly improved test accuracy and a slightly decreased test AUC. However, it does not solve the potential problem of overfitting as train performance is still significantly better than test performance. In fact, train performance on the reduced random forest is better than the random forest with a full set of variables - the reduced complexity of the model enabled it to have a lower bias on the training set.
XGBoost
For completeness, we also consider the xgboost library for gradient boosted decision trees. Gradient Boosted Decision Trees. The XGBoost package, introduced in (Chen, Tianqi and Guestrin, Carlos 2016) is a variant of Gradient Boosted Decision Trees (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome 2009, 353–74)
dmy_train <- dummyVars("~.", data = df_train[,-length(df_train)])
dmy_test <- dummyVars("~.", data = df_test[,-length(df_test)])
X_train <- as.matrix(data.frame(predict(dmy_train,df_train)))
X_test <- as.matrix(data.frame(predict(dmy_test,df_test)))
y_train = as.integer(as.matrix(df_train$target))
y_test = as.integer(as.matrix(df_test$target))
bst <- xgboost(data = X_train, label=y_train, max_depth = 2, nround = 10,
verbose=0,
objective='binary:logistic',
eval_metric="logloss")
pred_train <- predict(bst, X_train, type="response")
pred_test <- predict(bst, X_test, type="response")
results["xgb",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["xgb",]
importance_matrix <- xgb.importance(model=bst)
xgb.plot.importance(importance_matrix)
# shapley values
xgboost::xgb.ggplot.shap.summary(X_test, model = bst,
target_class = 1, top_n = 20) # Summary plot
CatBoost
We also consider briefly examine the catboost library for gradient boosted decision trees, given its popularity on machine learning competitions such as Kaggle. The CatBoost library has the advantage of learning a target encoding for categorical variables. From an implementation point of view, this may reduce the preprocessing that may be required.For details on CatBoost, refer to (Prokhorenkova, Liudmila and Gusev, Gleb and Vorobev, Aleksandr and Dorogush, Anna Veronika and Gulin, Andrey 2017)
X_train = df_train[,c(cat_feats, cont_feats)]
y_train = as.integer(df_train$target)
X_test = df_test[,c(cat_feats, cont_feats)]
y_test = as.integer(df_test$target)
pool <- catboost.load_pool(X_train, y_train, cat_features = cat_feats)
model <- catboost.train(pool, params=list(depth = 8, iterations = 10,
loss_function='Logloss', verbose=0))
pred_train <- catboost.predict(model, catboost.load_pool(X_train), prediction_type = 'Probability')
pred_test <- catboost.predict(model, catboost.load_pool(X_test), prediction_type = 'Probability')
# feature importance
feat_importance <- catboost.get_feature_importance(model, pool)
importances <- data.frame(feat_importance[order(feat_importance, decreasing=FALSE),])
importances$features = rownames(importances)
names(importances) <- c("importance","features")
importances$features <- factor(importances$features, level=importances$features)
ggplot(importances, aes(x=importance, y=features)) +
geom_bar(stat="identity") +
ggtitle("CatBoost Feature Importance")
#Shapley Values
data_shap_tree <- catboost.get_feature_importance(model, pool = pool,
type = "ShapValues")
data_shap_tree <- data.frame(data_shap_tree[, -ncol(data_shap_tree)])
names(data_shap_tree) = names(df[, c(cat_feats, cont_feats)])
ggplot(stack(data_shap_tree), aes(x = ind, y = values)) +
geom_point(aes(color = values)) + coord_flip() +
ggtitle("Shapley Values by variable") + scale_color_viridis_c()
results["catboost",] <- diagnosis(pred_train, pred_test, df_train$target, df_test$target)
results["catboost",]Results and Evaluation
results
results[order(-results$test_auc, -results$test_acc),]Although the tree-based methods are less interpretable and more “black box” to some extent, we can make use of model interpretability techniques.
Bibliography
Appendix
set.seed(2021)
cv_splits <- rsample::vfold_cv(df_train, strata = target, v= 3)
mod <- logistic_reg(penalty = tune(),
mixture = tune()) %>%
set_engine("glmnet")
glmnet_recipe <- recipe(target~.,data = df_train) %>%
step_dummy(all_nominal(), -all_outcomes())
glmnet_workflow<- workflow() %>%
add_recipe(glmnet_recipe) %>%
add_model(mod)
glmn_set <- parameters(penalty(range = c(-5,1), trans = log10_trans()),
mixture())
glmn_grid <-
grid_regular(glmn_set, levels = c(7, 5))
ctrl <- control_grid(save_pred = TRUE, verbose = TRUE)
glmn_tune <-
tune_grid(glmnet_workflow,
resamples = cv_splits,
grid = glmn_grid,
metrics = metric_set(roc_auc),
control = ctrl)
best_glmn <- select_best(glmn_tune, metric = "roc_auc")
glmnet_model <-
logistic_reg() %>%
set_engine("glmnet", seed=2021) %>%
set_mode("classification")
glm3 <- fit(glmnet_workflow, data = df_train)xgb_spec <- boost_tree(
trees = 100,
tree_depth = tune(),
learn_rate = tune() ## step size
) %>%
set_engine("xgboost") %>%
set_mode("classification")
xgb_grid <- grid_latin_hypercube(
tree_depth(),
learn_rate(),
size = 5
)
xgb_wf <- workflow() %>%
add_formula(target ~ .) %>%
add_model(xgb_spec)
set.seed(123)
vb_folds <- vfold_cv(df_train, strata = target)
doParallel::registerDoParallel()
set.seed(234)
xgb_res <- tune_grid(
xgb_wf,
resamples = vb_folds,
grid = xgb_grid,
control = control_grid(save_pred = TRUE)
)
best_auc <- select_best(xgb_res, "roc_auc")
final_xgb <- finalize_workflow(
xgb_wf,
best_auc
)
library(vip)
final_xgb %>%
fit(data = vb_train) %>%
pull_workflow_fit() %>%
vip(geom = "point")
final_res <- last_fit(final_xgb, df_split)
collect_metrics(final_res)xgb_res %>%
collect_metrics() %>%
filter(.metric == "roc_auc") %>%
select(mean, mtry:sample_size) %>%
pivot_longer(mtry:sample_size,
values_to = "value",
names_to = "parameter"
) %>%
ggplot(aes(value, mean, color = parameter)) +
geom_point(alpha = 0.8, show.legend = FALSE) +
facet_wrap(~parameter, scales = "free_x") +
labs(x = NULL, y = "AUC")"https://curso-r.github.io/treesnip/index.html"
library(treesnip)